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Abstract 

In random matrix theory, the spacing distribution functions p( n \s) 
are well fitted by the Wigner surmise and its generalizations. In this 
approximation the spacing functions are completely described by the 
behavior of the exact functions in the limits s — > and s — > oo. 
Most non equilibrium systems do not have analytical solutions for 
the spacing distribution and correlation functions. Because of that, 
we explore the possibility to use the Wigner surmise approximation 
in these systems. We found that this approximation provides a first 
approach to the statistical behavior of complex systems, in particular 
we use it to find an analytical approximation to the nearest neighbor 
distribution of the annihilation random walk. 

Keywords: Systems out of equilibrium, random matrices, Wigner sur- 
mise. 



1 Introduction 

In random matrix theory the analytic expressions for the spacing distribu- 
tion functions of eigenvalues in the circular and Gaussian orthogonal 
ensembles (COE and GOE respectively) in the limit of large matrices are 
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given in terms of the eigenvalues and eigenf unctions fi(x) of the following 
integral equation, see Ref. pQ: 



Hfi(x) = / e^^mdy. (1) 



i 



The spacing distributions are calculated explicitly by using 
E(2r, s) = H(l " A 2l ) Yl II ( >< [1 " + " " " + b jr )] , 

i=0 0<ji <j 2 <•••<> i=l ^ 

(2) 

S(2r-l, a ) = II(l-A*) £ ntl^F ) x ^ + --- + & >)' ( 3 ) 



where 
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h i = hjiX) J f2j(x)dx/ j ^ f2j(x)dx, (4) 
A, = s|^| 2 /4, (5) 



and 



j2 n 

= ^5>- i + «*)■ ( 6 ) 

S i=o 

These expressions are difficult to manage, however in Ref. [2j, the authors 
find an excellent approximation for spacing distributions p( n '(s) from their 
well-known behavior in the limits s — > and s — > 00. This approximation 
is easy to use and provide an excellent fit to the exact distributions. We 
will use this approximation many times in this paper, because of that, we 
summarize now its most important aspects. 

By definition, p( n \s) is the probability density that an interval of length 
s which starts at a level contains exactly n levels and the next, the n + 1 
level, is in [s, s + ds). In the same way, let F' n '(s) be the probability that an 
interval of length s which starts at a level, contains n levels. By using this 
definition we can write 
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(p^^-p^is'^ds'. 



(7) 



Additionally, let r^ n \s) be the probability density that an interval [0, s] which 
starts at a level at s = is limited by a level on its right side, under the 
condition that there are exactly n levels in the interval (0,s), i.e., r^ n \s) is 
the conditional probability 



r (n) (s) 



P 



in) I 



(8) 



this probability is called level repulsion function. Following Ref. [2], in the 
limit s — > 0, this equation can be written as 



p( n \s) 



r in) (s) 



P 



in-l)(J 



s')ds'. 



(9) 



In the GOE ensemble the matrix elements are chosen using a Gaussian 
distribution, this fact suggest that p( n \s) decays as Gaussian function. The 
appropriate function for fit is 



p {n \s) 



A „a n -B n s 2 

si n h e , 



(10) 



under the surmise r^ n \s) 
» 



with s — > 0. Additionally, the functions 



p^ n '(s) satisfy the normalization conditions 



POO 

/ p {n) (s)ds = l, 
Jo 



and 



sp {n) (s)ds = 1. 



(11) 
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By using the surmise for the level repulsion and the normalization condi- 
tions, is straightforward to find [2] 



R (a„+l)/2 
A - 0±^ 

r(f) : 
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r(f + i) 
+ i)r (atfcl) 



(13) 
(14) 
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where 



a. 



n + 



(n + l)(n + 2) 



(15) 



Then, the approximate spacing distribution functions p( n '(s) are given 
explicitly by 



(n+l) 



r(s^ti) an - 



r(gf+i) 
( „+i )r (£**+!)" 



(16) 



The result obtained for a n coincides with the results obtained by using 
the exact expression for the spacing distribution functions, see Ref. p]. No- 
tice that the approximate spacing distributions functions are characterized 
by the level repulsion, normalization condition, scaling condition for the av- 
erage spacing and Gaussian decay. This approximation is called generalized 
Wigner surmise and provides a very good approximation for p( n '(s), because 
it reproduce not only the distributions behavior in the limits s — > and 
s — > oo, but also reproduce their global behavior, as we can see in figure [TJ 
In particular the function with n = is called Wigner distribution. This 
fit allow us calculate also the approximate pair correlation distribution g(r). 
For this purpose we use 



g(r) = J2p {n \r) 



(17) 



n=0 



then 



9{r) 
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(n + l) 
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(18) 



In figured] we can see that this is a good approximation for g(r), however, it 
is not as useful as the Wigner surmise for p(°> (s) because the exact expression 
for g(r) is well known and easy to use, see Ref. pp. 

In Ref. |3J the authors study the statistical behavior of several out of equi- 
librium domain systems which evolve with formation of domains which grow 
in time. For intermediate times where the size of the domains is much smaller 
than the total size L of the system, the domain size distribution exhibit a 
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dynamic scaling. The authors studied the statistical properties of these do- 
mains in the scaling regime. They found that the statistical behavior of those 
is similar to the one in random matrices, for example, the nearest neighbor 
distribution p^(s) of several out of equilibrium domain systems is well fitted 
by the Wigner surmise which also describe closely the distribution p(°)(s) in 
the case of the circular and Gaussian orthogonal ensembles in random ma- 
trix theory (actually this distribution is exact in the case of 2 x 2 matrices). 
However, the next distributions (n > 0) for these systems are different from 
their counterpart in random matrix theory. Another important aspect is the 
pair correlation function g(r) which, in COE and GOE ensembles and the 
coalescing random walk and interacting random walk does not have any os- 
cillation but in other systems g(r) describe one oscillation near to r = 1. For 
more information see Refs. [3], HI [5j [6] . In most of the non equilibrium domain 
systems, the main problem is the absence of analytical expressions for the 
spacing and correlation functions. Then, the question is: can the generalized 
Wigner surmise provide a good approximation for p^^s) and g(r) in the 
domain systems as it happens with the random matrix ensembles? 




s , r 



Figure 1: Comparison between the generalized Wigner surmise and the COE 
ensemble. In the simulation we took 20000 matrices of size 200 x 200. 
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2 Wigner surmise for domains systems 



For all systems considered in this paper (s) is well described by the Wigner 
distribution, because of that and following the method used in the random 
matrix theory we propose the next model 

f 1 for n = , , 

an \ h(n) for n > 1 1 ' 

with h(n) is a function to determine. The spacing distribution functions in 
this model are given by 

P An&e-*** ifn>l ' 



using (fTTj) and ffT2l) . we find 



0n 



and 



p / 2+On \ 

^ = '(l + n)?(^l ' (22) 



2.1 Independent interval approximation model (HA) 

The independent intervals are used as an approximate solution in many equi- 
librium and non equilibrium systems [3j [8] in order to find analytical re- 
sults. In this approximation, p(")(s) is given by the convolution product of 
n + 1 nearest neighbor distribution factors, because of that, the spacing dis- 
tribution functions can be calculated by using the Laplace transformation, 
see Ref. [7]. In particular, in Ref. [3] the IIA is used to find an approximate 
model for the statistical behavior of two non equilibrium systems which will 
be explained in next sections. 



2.1.1 Independent interval model for small values of s 

In Ref. [3] the authors choose p^(s) equal to the Wigner distribution. In 
order to apply the method of the last section, we need to know the behavior 
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of p( n \s) for small and large values of s. For the first region we expand the 
Wigner distribution in power series 

p(o)( a) = ^ e -i- a = | a (i-p + ...), (23) 

then, to the first order, the nearest neighbor distribution p(°\s) has a lineal 
behavior, given by 

p (0) (s)ocs. (24) 

In the same limit s — > 0, by using the independent interval approximation 
for arbitrary values of n, we have 

J 0<xi<X2---<x„<s 

which can be evaluated by using the Laplace transform 

^(t) « 7^1) > (26) 

and then, taking its inverse 

p {n) (s) ocs 2n+1 . (27) 

As consequence, in the IIA case the exponent a n depends linearly on n 

a n = 2n + l. (28) 

By using this result it is possible to determine the behavior of the level 
repulsion function for s — > 0. Following Ref. [2] we have 

r (n) (s) oc s /(n) , (29) 

where f(n) is the function to determine. By using equation we can write 



j9 W (s)ocs /w / p^-^(s')^, (30) 

then 



o 



p^( S )oc^ + - +/ ^ +n . (31) 
By comparing (|27|) with (!3T|) is straightforward to find 

/(n) = 1, (32) 

for all n > 0, as a consequence 

r (n) (s) ocs, s^O, (33) 
then, the level repulsion does not depend on n as it happens in the COE/GOE 



case. 
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2.1.2 Independent interval model for large values of s 

Now, we need the behavior of p^(s) for large values of s. The exact expres- 
sion for p( n \s) is 

p (n) (s) = f p i0 \x 1 )p i0 \x2-x 1 )---p {0 \s-x n )dx 1 ---dx n . (34) 

J 0<xi<X2---<x n <s 

In our case p(°'(s) is given by the Wigner surmise, then 
P (n \s) = Hpj / xie~^ Xl ' 2 ■ ■ ■ (s - x n )e~^ {s ~ Xn)2 dx 1 ■ ■ ■ dx n . 

Jo<xi<x 2 -<x n <s 

(35) 

We can calculate the behavior of these functions for arbitrary values of n 
in this limit s — > oo as we show next. From Ref. [3] we know that at least 
the first two spacing distribution functions decay like Gaussian functions, 
then, we assume that for arbitrary values of n these functions have the form 
Pasy{s) = M n s' yn e~ Nns2 in the limit s — > oo. In order to eliminate the integrals 
in equation (135]) we use the Laplace transformation 

ft n \l) = (l - Ze /2/w erfc ^ + , (36) 

where erfc(z) = (2/ y^r) e~* 2 dt is the complementary Gaussian error func- 
tion. In the same way we take the Laplace transform in p^asy- Additionally, 
we expand both transformations in Taylor series around I = 0. Let be Zj 



the j coefficient in the expansion of equation (|36jl and Yj is the one for 
the Laplace transform of Pasy(s). We find that the coefficients of both ex- 
pansions satisfy the relation Yi/Zi = Yj/Zj in the limit i, j — > oo. By using 
this method we can find M n , N n and (3 n . If fact we find that N n = j- and 
In — n + 1. In general, if we know p^(s) we can calculate the asymptotic 
behavior of p( n \s) under the assumption that the IIA is valid for s — ► oo, 
but, as we will see in next sections, this is not true always. 

In the figure [2] we compare the exact statistical behavior of IIA with the 
generalized Wigner surmise, i.e., with a fit developed by using the behavior 
of p^ (s) in the limits s — » and s — > oo, because of that from now on we will 
call it local fit. Also, we compare the global fit which was developed by using 
equations f|T9|) to (1221) and the complete behavior of p( n \s) in the interval 
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[0, oo]. By using the values of a n found in the global fit, we developed a new- 
fit to determine the global behavior of a n , explicitly in this case we have 



this result is close to the exact exponent fl28|) . even when we use wrong 
functions in the fit; for example, the exact result for pW(s) is, see Ref. [3] 



which is very different form our surmise, however, both functions (120]) and 
( I38p have the same type of behavior in the limits s — > and s — > oo. Equation 
( I18p for the correlation function it is still valid in both cases, global and local 
fit, we only must use equation ( 128]) and ( 1371) respectively. The main problem 
in the global fit approximation it is the use of not integer exponents in the 
level repulsion. Figure [3] show the differences between the three cases for 
small values of s, naturally in this region the graph of the global fit is not 
parallel to graph of the exact result as it actually happens in the local fit 
approximation. In figure H] we can see the linear behavior of p^ n '(s) in limit 
s — > oo, which implies that the distribution functions decay like a Gaussian 
function as it was to be expected. 

2.2 Coalescing random walk (CRW) 

In the coalescing random walk the particles describe independent random 
walks along a one dimensional lattice and they are subjected to the reaction 
A + A — > A. This system is well studied [U EJ EH [TT] and its analytical 
solution is well know, because of that is used as approximation to more 
complex systems. Let q( n '(s) be the conditional probability that given one 
particle its next neighbor is at a distance of s. From its definition q( n \s) is 
given by 



a n = 1.826871 + 0.9954, 



(37) 



V 





(39) 



with 





gn E (n V( Xu y u ... , X w _i, ?7 n _i) 

dx 1 ■ ■■dx n _idy n _ l 



yi=X2,-- ,y n -i=x. 



, (40) 
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Figure 2: Comparison between the exact statistical behavior of IIA, the 
generalized Wigner surmise (local fit) and the global fit. 




Figure 3: Log-Log graphic for the spacing distribution functions for IIA. 



10 




Figure 4: Asymptotic behavior for s — > oo of the spacing distribution func- 
tions. 



(271-1)!! 

E {n \x 1 ,y 1 , ■ ■ ■ ,x n ,y n ,t) = ^ & P E {1) ( z i, P , Z 2, P , t) ■ ■ ■ E {1 \z 2n -i, P , z 2n , P , t), 

P =i 

(41) 

where Zx tP , z 2tP , z 2niP symbolize an ordered permutation, p, of the variables 
xi,y h ...,x n ,y n , such that 

Z 1, P < Z 2, P , Z 3, P < ^4,p, - - - , Z2n-l,p < Z 2n ,p, (42) 

and 

z l, P < Z 3, P < Z 5, P ■ ■ ■ < Z 2n-l,p- 

(43) 

The function E^(xi, y±, t) is the probability that from x\ to y\ the lattice 
is empty at time t. Then it is possible generate the complete solution for 
the CRW from E^\xi, yi,t), which is given by the solution of the diffusion 
equation under the suitable boundary conditions (see Ref. jl]). In fact, the 
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exact expression for this function is 

EW(x uyi ,t)=eric 



Vi - xi 



V8Di 



(44) 



with D the diffusion constant and t the time, for additionally information 
see Ref. [1]. For practical purposes, the solution given by equations (1591) to 
f !44p is hard to evaluate for arbitrary values of n but it can be evaluated in 
the limit s — > using Taylor series. The case n = is trivial, the Taylor 
expansion for equation (Jill is 



E^(x u y u t) = l- 
then 



yi-xi , {yi-xi) 



2 7 r( J Dt) 1 /2 24 v / 27r( J Dt) 3 / 2 640V27r( J Dt) 5 /2 
d 2 



q(°\x 2 ,xi) = lu {2 \x 1 ,x 2 ) ■■ 

(0)/ \ x 2 — X\ 



dxidyi 

{x 2 - ^i) 3 

4 v / 2¥(L>t) 3 / 2 32 v / 2T( J Dt) 5 / 2 



+0(x,y) 7 , 
(45) 

(46) 



yi=X2 

+ 0(xf. 



(47) 



Making the variable change s 



V2TTDt 



and taking into account that p(°>(s) 



2nDt q(°>(x 2 , x%), the product Dt disappears (dynamical scaling) in the above 
equation. Then, to first order, we have 



P 



(0), 



STY 



+ 0(s) 



(48) 



For small values of s, p(°\s) has a linear behavior, i.e., «o = 1. The case 
n = 1 is more complicated, in fact we have 



uj {3) (xt, x 2 , x 3 ) = - 



d 3 



dx\dx 2 dy> 



■EW(x 1 ,y 1 ,x 2 ,y2,t) 



(49) 



5/1=^2, J/2 =Z'3 



where 



E {2) (x 1 ,y 1 ,x 2l y 2 ,t) = E (xi,y±, i) E (x 2 , y 2 , t) 

+ E(xi,y 2 ,t) E(y 1 ,x 2 ,t) 
- E(x 1 ,x 2 ,t) E( yi ,y 2 ,t) 



(50) 
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then 

(3)/ \ ( X 2 — Xx)(x3 — Xi)(x3 — X 2 ) , , 4 , . 

uW(x u X2,x 3 ) = S2n(Dt) 3 + °( x ) > ( 51 ) 



in that way q^(x\, x$, t) is given by 

ff (1) (*s,*i,f) = r {X2 ~ Xl) t 3 7"i {X3 ~ X2) d^ + Q(xf. (52) 



Xl 



327r(Dt) 3 



Integrating 



*' = fe^ + °W' (53) 

Using again the variable change, it is straigthfoward to find 

p a)( a ) = ^ + 0(x) 6 , (54) 

we conclude that a± = 4. In general for an arbitrary value of n, we find that 
the first term in the expansion is 

fin fX3 ^ ^ 

g (n) (xi,x n ,t) oc / ■■■ / ] [ (xj - Xi)dx 2 ■ ■ ■ dx n _ u (55) 

therefore, the above equation has (n+ l)(n + 2)/2 different factors which im- 
plies that the integrand is proportional to x\ n+1 ^ n+2 ^ 2 . Making the integral 
and the usual variable change, the final expression for small values of s is 
proportional to s ( n + 1 )( n + 2 )/ 2 + ri ) explicitly, we have 

(n + l)(n + 2) 
a n = nA . (56) 

This is the same result reported in Ref. [2] for the GOE/COE case and 
coincides with the partial result presented in Ref. pj] for the CRW. We made 
again both fits, global and local. The global fit was made with the data from 
our simulation where we use a lattice with 1000 sites and 500 particles in 
t = 0. The data to build the histograms was taken at three different times 
T = 50, T = 100 and T = 200 over 50000 realizations. In this case the global 
fit is not as accurate as in the IIA case as we can see in figure but it still is 
a good approximation. We use again equations ffl9l to (122]) ; and additionally 
we supposed a Gaussian decay (/3 = 2). The global fit gives 
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a n = 2.8688n + 0.8621. (57) 

The global fit gives an erroneous exponent which depend linearly with n, 
this result it does not coincide with the analytical result (|56l) . where, a n is 
a quadratic function of n. The local fit it is very different from the simu- 
lation results and coincides with the statistical behavior of the COE/GOE 
ensembles. 
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Figure 5: Comparison between the statistical behavior of CRW, the global 
fit and the local fit. 

Although the global and local fit models are approximate, we can use 
them as a good approximations in some cases. For example, in Ref. [TT] the 
authors find an exact relation for the nearest neighbor distribution pmn(s) 
in the annihilation random walk in terms of p( n )(s) of the coalescing random 
walk. Explicitly, they found 

n>0 

In order to test the validity of our approximations, we implement a simu- 
lation for the annihilation random walk for a one dimensional lattice with 
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Figure 7: Statistical behavior of CRW for small values of s. 
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2000 sites, 100 particles at t = over 20000 realizations, the histogram was 
build by using three times T = 1000, T = 1500 and T = 2000. By using 
the global and the local fit for the distribution functions p^(s) of the CRW, 
with equation (158]) . we find two analytical models for the annihilation ran- 
dom walk. We can see in figure [8] that the global and local fits provides a 
good approximation for the nearest neighbor distribution of the annihilation 
random walk. Additionally, figure M compare global and local fit with the 
asymptotic result pann(s) ~ 1.8167e~ L3062s given in Ref. [TT] , 



0,0 0,5 1.0 1.5 2,0 2,5 3,0 
■ > i i ■ - * - , 




s 



Figure 8: Approximation for pann(s) by using global and local fits. 



2.3 Spin System 

This system was introduced in Ref. [5] , where the authors consider a chain of 
L Ising spins with nearest neighbor ferromagnetic interaction J. The chain is 
subject to spin-exchange dynamics with a driving force E that favors motion 
of up spins to the right over motion to the left. In this case we do not 
have an analytical solution for the spacing distribution functions, because of 
that, we must start exploring numerically the behavior of p( n \s) for small 
and large values of s. In figure [TU1 we can see the linear behavior of the 
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Figure 9: Asymptotic behavior of p a nn(s). 



spacing distribution function for s — > 0. Using values in this region we 
develop a fit which suggest that a% — 3 and a<i = 6 approximately. Naturally 
«o = 1, however it is very difficult to know using this method the next 
exponents because it is not possible develop a numerical simulation with 
enough precision. 

Curiously, these exponents for n = 0, 1, 2 are given by the equation 

a „= ( " + 1) 2 ( " + 2) , (59) 

which is very similar to its counterpart in COE and GOE cases. For s — > oo, 
p( n \s) decay like a Gaussian function as we can see in figure fTTl In this case 
the global fit gives 

a n = 1.270n + 0.920. (60) 

In figure [12], we show the results given by equations (jl9~l) to (122"]) for the 
global fit in comparison with the simulation results which was made with a 
lattice with 1000 sites, equal number of spins up and down taken at two times 
t = 34 and t = 48 to build the histograms. The result for g(r) is very good 
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Figure 10: Log-Log graphic for the spacing distribution functions of spin 
system. 




Figure 11: Asymptotic behavior of the spin system. 
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with a maximum error of 2.5%. Unfortunately this approximation is not 
good enough for p( n >(s) but at least it reproduce qualitatively the behavior 
of the real functions for s — > oo. The local fit gives terrible results as it 
happens in the CRW case. 



Figure 12: Comparison between the statistical behavior of the spin system 
and the global fit. 



2.4 Gas System 



This system was originally studied in [Bj. There, the authors studied the 
biased diffusion of two species in a fully periodic 2 x L rectangular lattice 
half filled with two equal number of two types of particles (labeled by their 
charge + or — ). An infinite external field drives the two species in opposite 
directions along the x axis (long axis). The only interaction between particles 
is an excluded volume constraint, i.e., each lattice site can be occupied at 
most by only one particle. As it happens in the spin system, we do not know 
an analytic solution for the spacing and pair correlation functions. We follow 
the same method used in the spin system. In figure [T3J we can see the linear 
behavior of p( n \s) which, by fit, give us a x = 3 and a 2 = 5 approximately, 
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and of course ctQ = 1. This fact suggest a linear behavior for a n given by 

a n = 2n + 1 (61) 

but again we could not find the next exponents with enough precision in 
order to validate the above equation. For s — ► oo we found that p^°\s) 
decays like a Gaussian function {(3 = 2), but for n > 0, we found that 
(3 is an indeterminate function of n. For example in figure [TH we can see 
the asymptotic behavior for two consecutive spacing distribution functions, 
the figure sugg est p = 2 for p(°)(s) and (3 ^ 2 for pW(s) as it happens in 
Ref. [12]. Because it is difficult determine the exact value of f3 from the 
graphics, we implement a linear regression to find which value of j3 give us a 
better "straight" line. With this method we find for example, that (3 = 2.6 
for n = 1, (3 — 3 for n = 5 and /3 = 3.2 for n = 8. In the linear regressions 
we took values between 5.5 > s > 2.5, 11.7 > s > 7 and 15.5 > s > 10 
respectively. Because of that, for the gas system we propose a model where 
j3 depends on n. In particular we choose (3 = 2.6 + 0.1(n — 1). With this 
model, the global fit gives 

a n = 1.016n + 0.788. (62) 

The results of the global fit are show in figure [151 again we find good fit 
for g(r) with a maximum error of 2% approximately but the agreement for 
p^(s) is not so good. Additionally, we include the first spacing distribution 
obtained with the local fit and our model for (3 n . 

3 Conclusion 

In COE and GOE ensembles, the spacing distribution functions p( n '(s) can be 
well described by using their behavior for small and large values of s (local 
fit) as it happens in IIA case, however, this is not true for more complex 
systems like CRW, spin and gas systems. This result was to be expected 
because in general the spacing distribution functions are characterized also 
by their inter medium behavior. In general, the global fit gives better results 
in comparison with the local fit but it fails to reproduce the level repulsion, 
in fact, gives non integer exponents. The level repulsion for the CRW has 
the same behavior that the circular and Gaussian orthogonal ensembles, i.e., 
both systems are equivalents for s — > 0. The numerical results suggest that 
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Figure 13: Log-Log graphic for the spacing distribution functions of gas 
system. 

the IIA and the gas system are also equivalents in that region. We find 
numerical evidence that the spacing distributions functions for gas system is 
described by a non universal function, in fact, they decay as M n s ln e~ Nn s>3n 
for n > 0, with f3 n an indeterminate function of n. In general the global and 
local fit provides a first approximation for p( n \s) and g(r), which can be used 
as a good approximation as it happens in the annihilation random walk case. 
These approximations also serve to classify the spacing distribution functions 
according to their level of repulsion and their decay functional form. 
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Figure 15: Comparison between the statistical behavior of the gas system 
and the global fit. 
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